Skyler Lewis 2025-01-02
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(habistat)
library(patchwork)
theme_set(theme_minimal())
See watershed_delineation.R for dependencies
make_flow_xw <- function(group_var) {
habistat::flowline_attr |>
select({{group_var}},
comid = comid,
# da_reach = da_area_sq_km_tot,
# 20241223: Switched to divergence-routed due to unexpected large vlaues in da_area_sq_km_tot, e.g. see comid 1889628
da_reach = da_area_sq_km_div,
pc_reach = da_ppt_mean_mm) |>
drop_na({{group_var}}) |>
group_by({{group_var}}) |>
mutate(outlet_comid = comid[which.max(da_reach)]) |>
mutate(da_outlet = da_reach[which(comid == outlet_comid)],
pc_outlet = pc_reach[which(comid == outlet_comid)]) |>
mutate(multiplier = (da_reach / da_outlet) * (pc_reach / pc_outlet)) |>
arrange({{group_var}}, -multiplier)
}
cv_watersheds_flow_xw <- make_flow_xw(watershed_level_3)
cv_watersheds_flow_xw |> usethis::use_data(overwrite=T)
## ✔ Setting active project to
## "C:/Users/skylerlewis/Github/swc-habitat-suitability-parallel".
## ✔ Saving "cv_watersheds_flow_xw" to "data/cv_watersheds_flow_xw.rda".
## ☐ Document your data (see <https://r-pkgs.org/data.html>).
cv_mainstems_flow_xw <- make_flow_xw(river_cvpia)
cv_mainstems_flow_xw |> usethis::use_data(overwrite=T)
## ✔ Saving "cv_mainstems_flow_xw" to "data/cv_mainstems_flow_xw.rda".
## ☐ Document your data (see <https://r-pkgs.org/data.html>).
cv_mainstem_groups_flow_xw <- make_flow_xw(river_group)
cv_mainstem_groups_flow_xw |> usethis::use_data(overwrite=T)
## ✔ Saving "cv_mainstem_groups_flow_xw" to "data/cv_mainstem_groups_flow_xw.rda".
## ☐ Document your data (see <https://r-pkgs.org/data.html>).
habistat::flowline_geom_proj |>
inner_join(cv_watersheds_flow_xw, by=join_by(comid)) |>
filter(comid %in% habistat::wua_predicted$comid) |>
ggplot() +
geom_sf(aes(color = multiplier)) +
geom_sf(data = habistat::flowline_geom_proj |>
filter(comid %in% cv_watersheds_flow_xw$outlet_comid) |>
st_line_sample(sample=1),
aes(color = 1)) +
scale_color_viridis_c(direction = -1) +
ggtitle("Flow Multipliers by Watershed")
habistat::flowline_geom_proj |>
inner_join(cv_mainstems_flow_xw, by=join_by(comid)) |>
filter(comid %in% habistat::wua_predicted$comid) |>
ggplot() +
geom_sf(aes(color = multiplier)) +
geom_sf(data = habistat::flowline_geom_proj |>
filter(comid %in% cv_mainstems_flow_xw$outlet_comid) |>
st_line_sample(sample=1),
aes(color = 1)) +
scale_color_viridis_c(direction = -1) +
ggtitle("Flow Multipliers by Mainstem")
habistat::flowline_geom_proj |>
inner_join(cv_mainstem_groups_flow_xw, by=join_by(comid)) |>
filter(comid %in% habistat::wua_predicted$comid) |>
ggplot() +
geom_sf(aes(color = multiplier)) +
geom_sf(data = habistat::flowline_geom_proj |>
filter(comid %in% cv_mainstem_groups_flow_xw$outlet_comid) |>
st_line_sample(sample=1),
aes(color = 1)) +
scale_color_viridis_c(direction = -1) +
ggtitle("Flow Multipliers by Mainstem Group")
# habistat::wua_predicted <-
# readRDS(here::here("data-raw", "results", "habistat::wua_predicted.Rds"))
# for backwards compatibility because habistat::wua_predicted is now pivoted
# load(here::here("data", "wua_predicted.Rda")) # this contains the ensemble predictions
wua_predicted_cv_watersheds <-
habistat::wua_predicted |>
inner_join(cv_watersheds_flow_xw,
by=join_by(watershed_level_3, comid)) |>
expand_grid(scaled = c(FALSE, TRUE)) |>
mutate(flow_cfs = if_else(scaled, flow_cfs * multiplier, flow_cfs)) |>
arrange(habitat, comid, flow_cfs) |> # , model_name,
group_by(habitat, comid) |> # , model_name,
mutate(across(c(wua_per_lf_pred, wua_per_lf_pred_SD, wua_per_lf_pred_SN),
\(v) if_else(scaled, NA, v))) |>
mutate(across(c(wua_per_lf_pred, wua_per_lf_pred_SD, wua_per_lf_pred_SN),
\(v) zoo::na.approx(v, x = flow_cfs, na.rm=F, rule=2))) |>
filter(scaled) |>
select(-scaled) |>
ungroup() |>
group_by(habitat, watershed_level_3, flow_idx) |> # , model_name,
summarize(wua_per_lf_pred = sum(wua_per_lf_pred * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_per_lf_pred_SD = sum(wua_per_lf_pred_SD * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_per_lf_pred_SN = sum(wua_per_lf_pred_SN * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_acres_pred = sum(wua_per_lf_pred * reach_length_ft, na.rm=T) / 43560,
wua_acres_pred_SD = sum(wua_per_lf_pred_SD * reach_length_ft, na.rm=T) / 43560,
wua_acres_pred_SN = sum(wua_per_lf_pred_SN * reach_length_ft, na.rm=T) / 43560,
.groups="drop") |>
inner_join(habistat::wua_predicted |>
group_by(flow_idx) |>
summarize(flow_cfs = first(flow_cfs), .groups="drop"),
by = join_by(flow_idx))
# wua_predicted_cv_watersheds <- wua_predicted_cv_watersheds |>
# pivot_wider(names_from = model_name, values_from = c(wua_per_lf_pred, wua_acres_pred)) |>
# mutate(wua_per_lf_pred = (wua_per_lf_pred_SD + wua_per_lf_pred_SN) / 2,
# wua_acres_pred = (wua_acres_pred_SD + wua_acres_pred_SN) / 2)
wua_predicted_cv_watersheds |> usethis::use_data(overwrite=T)
## ✔ Saving "wua_predicted_cv_watersheds" to
## "data/wua_predicted_cv_watersheds.rda".
## ☐ Document your data (see <https://r-pkgs.org/data.html>).
wua_predicted_cv_mainstems <-
habistat::wua_predicted |>
inner_join(cv_mainstems_flow_xw,
by=join_by(river_cvpia, comid)) |>
expand_grid(scaled = c(FALSE, TRUE)) |>
mutate(flow_cfs = if_else(scaled, flow_cfs * multiplier, flow_cfs)) |>
arrange(habitat, comid, flow_cfs) |> # , model_name,
group_by(habitat, comid) |> # , model_name,
mutate(across(c(wua_per_lf_pred, wua_per_lf_pred_SD, wua_per_lf_pred_SN),
\(v) if_else(scaled, NA, v))) |>
mutate(across(c(wua_per_lf_pred, wua_per_lf_pred_SD, wua_per_lf_pred_SN),
\(v) zoo::na.approx(v, x = flow_cfs, na.rm=F, rule=2))) |>
filter(scaled) |>
select(-scaled) |>
ungroup() |>
group_by(habitat, river_cvpia, flow_idx) |> # , model_name,
summarize(wua_per_lf_pred = sum(wua_per_lf_pred * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_per_lf_pred_SD = sum(wua_per_lf_pred_SD * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_per_lf_pred_SN = sum(wua_per_lf_pred_SN * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_acres_pred = sum(wua_per_lf_pred * reach_length_ft, na.rm=T) / 43560,
wua_acres_pred_SD = sum(wua_per_lf_pred_SD * reach_length_ft, na.rm=T) / 43560,
wua_acres_pred_SN = sum(wua_per_lf_pred_SN * reach_length_ft, na.rm=T) / 43560,
.groups="drop") |>
inner_join(habistat::wua_predicted |>
group_by(flow_idx) |>
summarize(flow_cfs = first(flow_cfs), .groups="drop"),
by = join_by(flow_idx))
# wua_predicted_cv_mainstems <- wua_predicted_cv_mainstems |>
# pivot_wider(names_from = model_name, values_from = c(wua_per_lf_pred, wua_acres_pred)) |>
# mutate(wua_per_lf_pred = (wua_per_lf_pred_SD + wua_per_lf_pred_SN) / 2,
# wua_acres_pred = (wua_acres_pred_SD + wua_acres_pred_SN) / 2)
wua_predicted_cv_mainstems |> usethis::use_data(overwrite=T)
## ✔ Saving "wua_predicted_cv_mainstems" to "data/wua_predicted_cv_mainstems.rda".
## ☐ Document your data (see <https://r-pkgs.org/data.html>).
# these are the mainstems grouped for CVPIA model comparison
wua_predicted_cv_mainstems_grouped <-
habistat::wua_predicted |>
inner_join(cv_mainstem_groups_flow_xw,
by=join_by(river_group, comid)) |>
expand_grid(scaled = c(FALSE, TRUE)) |>
mutate(flow_cfs = if_else(scaled, flow_cfs * multiplier, flow_cfs)) |>
arrange(habitat, comid, flow_cfs) |> # , model_name,
group_by(habitat, comid) |> # , model_name,
mutate(across(c(wua_per_lf_pred, wua_per_lf_pred_SD, wua_per_lf_pred_SN),
\(v) if_else(scaled, NA, v))) |>
mutate(across(c(wua_per_lf_pred, wua_per_lf_pred_SD, wua_per_lf_pred_SN),
\(v) zoo::na.approx(v, x = flow_cfs, na.rm=F, rule=2))) |>
filter(scaled) |>
select(-scaled) |>
ungroup() |>
group_by(habitat, river_group, flow_idx) |> # , model_name,
summarize(wua_per_lf_pred = sum(wua_per_lf_pred * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_per_lf_pred_SD = sum(wua_per_lf_pred_SD * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_per_lf_pred_SN = sum(wua_per_lf_pred_SN * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_acres_pred = sum(wua_per_lf_pred * reach_length_ft, na.rm=T) / 43560,
wua_acres_pred_SD = sum(wua_per_lf_pred_SD * reach_length_ft, na.rm=T) / 43560,
wua_acres_pred_SN = sum(wua_per_lf_pred_SN * reach_length_ft, na.rm=T) / 43560,
.groups="drop") |>
inner_join(habistat::wua_predicted |>
group_by(flow_idx) |>
summarize(flow_cfs = first(flow_cfs), .groups="drop"),
by = join_by(flow_idx))
# wua_predicted_cv_mainstems_grouped <- wua_predicted_cv_mainstems_grouped |>
# pivot_wider(names_from = model_name, values_from = c(wua_per_lf_pred, wua_acres_pred)) |>
# mutate(wua_per_lf_pred = (wua_per_lf_pred_SD + wua_per_lf_pred_SN) / 2,
# wua_acres_pred = (wua_acres_pred_SD + wua_acres_pred_SN) / 2)
wua_predicted_cv_mainstems_grouped |> usethis::use_data(overwrite=T)
## ✔ Saving "wua_predicted_cv_mainstems_grouped" to
## "data/wua_predicted_cv_mainstems_grouped.rda".
## ☐ Document your data (see <https://r-pkgs.org/data.html>).
Plotting the output
drainage_areas <-
habistat::cv_watersheds |>
group_by(watershed_level_3) |>
summarize() |>
transmute(river_group = watershed_level_3,
area_ac = st_area(geometry) |>
units::set_units("acres") |>
units::drop_units()) |>
st_drop_geometry()
watershed_labels <-
drainage_areas |>
transmute(river_group,
label = str_glue("{river_group}\n{format(round(area_ac, -3), big.mark=',')} ac")) |>
deframe()
wua_predicted_cv_watersheds |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_level_3, scales="free") +
geom_ribbon(aes(ymin = wua_acres_pred_SD, ymax = wua_acres_pred_SN, fill = habitat), alpha=0.33) +
geom_line(aes(y = wua_acres_pred, color = habitat)) +
scale_x_log10() +
theme(legend.position = "top", panel.grid.minor = element_blank()) +
scale_color_brewer(palette = "Dark2", aesthetics = c("color", "fill")) +
ylab("Predicted Total Habitat (acres)") + xlab("Flow at Outlet (cfs)") +
ggtitle("Watersheds (incl. tributaries)")
wua_predicted_cv_mainstems |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~river_cvpia, scales="free") +
geom_ribbon(aes(ymin = wua_acres_pred_SD, ymax = wua_acres_pred_SN, fill = habitat), alpha=0.33) +
geom_line(aes(y = wua_acres_pred, color = habitat)) +
scale_x_log10() +
theme(legend.position = "top", panel.grid.minor = element_blank()) +
scale_color_brewer(palette = "Dark2", aesthetics = c("color", "fill")) +
ylab("Predicted Total Habitat (acres)") + xlab("Flow at Outlet (cfs)") +
ggtitle("Mainstems")
wua_predicted_cv_mainstems_grouped |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_ribbon(aes(ymin = wua_acres_pred_SD, ymax = wua_acres_pred_SN, fill = habitat), alpha=0.33) +
geom_line(aes(y = wua_acres_pred, color = habitat)) +
scale_x_log10() +
theme(legend.position = "top", panel.grid.minor = element_blank()) +
scale_color_brewer(palette = "Dark2", aesthetics = c("color", "fill")) +
ylab("Predicted Total Habitat (acres)") + xlab("Flow at Outlet (cfs)") +
ggtitle("Mainstems (grouped by watershed)")
Versus the “naive” method
wua_predicted |>
filter(!is.na(river_group)) |>
group_by(habitat, river_group, flow_idx, flow_cfs) |>
summarize(wua_per_lf_pred = sum(wua_per_lf_pred * reach_length_ft, na.rm=T) / sum(reach_length_ft, na.rm=T),
wua_acres_pred = sum(wua_per_lf_pred * reach_length_ft, na.rm=T) / 43560, .groups="drop") |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_line(aes(y = wua_acres_pred, color = habitat)) +
scale_x_log10() +
theme(legend.position = "top", panel.grid.minor = element_blank()) +
scale_color_brewer(palette = "Dark2", aesthetics = c("color", "fill")) +
ylab("Predicted Total Habitat (acres)") + xlab("Flow (cfs)") +
ggtitle("Example incorrect method of aggregating - do not use")
Moved this over from model_cleaned.Rmd
mainstems_comid <-
read_sf(file.path("/vsizip", here::here("data-raw", "source", "rearing_spatial_data", "nhdplusv2_comid_habitat_xw.shp.zip"))) |>
janitor::clean_names() |>
st_zm() |>
st_transform(st_crs(habistat::flowline_geom_proj)) |>
mutate(length_ft = st_length(geometry) |> units::set_units("ft") |> units::drop_units()) |>
filter(str_detect(habitat, "rearing")) |>
left_join(habistat::flowline_attr |>
select(comid, hqt_gradient_class, river_group), by=join_by(comid)) |>
filter(!(watershed %in% c("Sacramento River", "San Joaquin River")))
mainstems <-
mainstems_comid |>
group_by(watershed, river, hqt_gradient_class) |>
summarize()
## `summarise()` has grouped output by 'watershed', 'river'. You can override
## using the `.groups` argument.
mainstems_comid |>
ggplot() +
geom_sf(aes(group=river, color=hqt_gradient_class)) +
theme(legend.key.height = unit(12, "point")) +
ggtitle("DSMhabitat reaches")
# http://cvpia-habitat-docs-markdown.s3-website-us-west-2.amazonaws.com/watershed/Regional_Approximation.html
# These are the watersheds that use regional approximation for instream rearing habitat
regional_approx_groups <-
c("Antelope Creek", "Bear Creek", "Big Chico Creek",
"Elder Creek", "Mill Creek", "Paynes Creek", "Stony Creek", "Thomes Creek")
# These are the watersheds that use regional approximation for instream spawning habitat
regional_approx_groups_spawning <-
c("Antelope Creek", "Bear Creek", "Big Chico Creek", "Cow Creek",
"Elder Creek", "Mill Creek", "Paynes Creek", "Stony Creek", "Thomes Creek")
# These are the watersheds that use regional approximation for spawning and instream rearing in the san joaquin basin
regional_approx_groups_sanjoaquin <-
c("Cosumnes River")
# These are the watersheds that use scaled proxies for floodplain habitat
deer_creek_fp_proxy <-
c("Antelope Creek", "Bear Creek", "Cow Creek", "Mill Creek", "Paynes Creek")
deer_creek_partial_model_scaled <- c("Big Chico Creek")
cottonwood_creek_fp_proxy <- c("Stony Creek", "Thomes Creek")
tuolumne_river_fp_proxy <- c("Calaveras River")
fp_proxy_groups <- c(deer_creek_fp_proxy, cottonwood_creek_fp_proxy, tuolumne_river_fp_proxy)
# These are reaches where suitability is already applied so we shouldn't be using DSMhabitat::apply_suitability to scale again. List is from the help docs for DSMhabitat::apply_suitability
suitability_already_applied <- c("Antelope Creek", "Battle Creek", "Bear Creek", "Cow Creek", "Deer Creek", "Mill Creek", "Paynes Creek", "Sacramento River", "Sutter Bypass", "Yolo Bypass", "North Delta", "South Delta")
# pretty sure that all the scaled watersheds are in this too...
suitability_already_applied <- unique(c(suitability_already_applied, fp_proxy_groups))
#remotes::install_github("CVPIA-OSC/DSMhabitat")
watersheds <- mainstems |> pull(watershed) |> unique()
watershed_name <- tolower(gsub(pattern = "-| ", replacement = "_", x = watersheds))
watershed_rda_name <- paste(watershed_name, "floodplain", sep = "_")
dsm_habitat_floodplain <- map_df(watershed_rda_name, function(watershed) {
df <- as.data.frame(do.call(`::`, list(pkg = "DSMhabitat", name = watershed)))
}) |>
select(river_group = watershed,
flow_cfs,
ends_with("_floodplain_acres")) |>
pivot_longer(cols = ends_with("_floodplain_acres"),
names_transform = \(x) str_replace(x, "_floodplain_acres", ""),
names_to = "run",
values_to = "suitable_ac") |>
#values_to = "floodplain_acres") |>
mutate(run = run |> factor(levels = c("FR", "LFR", "WR", "SR", "ST"),
labels = c("fall", "late fall", "winter", "spring", "steelhead")),
hab = "floodplain" |> factor(levels = c("spawn", "fry", "juv", "adult", "floodplain"))) |>
mutate(river_group = if_else(river_group == "Consumnes", "Cosumnes River", river_group))
#floodplain_acres_suitable = if_else(river_group %in% suitability_already_applied,
# floodplain_acres,
# DSMhabitat::apply_suitability(floodplain_acres * 4046.86) / 4046.86))
dsm_habitat_instream <- map_df(paste(watershed_name, "instream", sep = "_"),
possibly(function(watershed) {
df <- as.data.frame(do.call(`::`, list(pkg = "DSMhabitat", name = watershed)))
}, otherwise = NULL)) |>
select(river_group = watershed,
flow_cfs,
ends_with("_wua")) |>
pivot_longer(cols = ends_with("_wua"),
names_transform = \(x) str_replace(x, "_wua", ""),
names_to = "run_hab",
values_to = "wua_per_1000ft") |>
separate_wider_delim(run_hab, names = c("run", "hab"), delim = "_") |>
mutate(wua_per_lf = wua_per_1000ft / 1000,
run = run |> factor(levels = c("FR", "LFR", "WR", "SR", "ST"),
labels = c("fall", "late fall", "winter", "spring", "steelhead")),
# TODO: use "labels" to align these with habistat names like "fall" instead of "FR"
hab = hab |> factor(levels = c("spawn", "fry", "juv", "adult", "floodplain"))) |>
drop_na() |>
select(-wua_per_1000ft)
dsm_habitat_instream_sqm <- map_df(paste(watershed_name, "instream", sep = "_"),
possibly(function(watershed) {
df <- as.data.frame(do.call(`::`, list(pkg = "DSMhabitat", name = watershed)))
}, otherwise = NULL)) |>
select(river_group = watershed,
flow_cfs,
ends_with("_sqm")) |>
pivot_longer(cols = ends_with("_sqm"),
names_transform = \(x) str_replace(x, "_sqm", ""),
names_to = "run_hab",
values_to = "suitable_sqm") |>
separate_wider_delim(run_hab, names = c("run", "hab"), delim = "_") |>
mutate(run = run |> factor(levels = c("FR", "LFR", "WR", "SR", "ST"),
labels = c("fall", "late fall", "winter", "spring", "steelhead")),
hab = hab |> factor(levels = c("spawn", "fry", "juv", "adult", "floodplain"))) |>
mutate(suitable_ac = suitable_sqm * (1 / 0.3048)^2 * (1 / 43560)) |>
select(-suitable_sqm)
# flow ranges
dsm_flows <- bind_rows(dsm_habitat_floodplain, dsm_habitat_instream) |>
group_by(river_group, flow_cfs) |>
summarize() |>
ungroup() |>
arrange(river_group, flow_cfs)
## `summarise()` has grouped output by 'river_group'. You can override using the
## `.groups` argument.
dsm_flow_ranges <-
dsm_flows |>
group_by(river_group) |>
summarize(min_flow_cfs = min(flow_cfs), max_flow_cfs = max(flow_cfs))
mainstems_comid |>
st_zm() |>
filter(comid %in% mainstems_comid$comid) |>
ggplot() +
geom_sf(aes(color=river_group)) +
theme(legend.key.height = unit(12, "point"))
cv_mainstems_spawning <-
habistat::cv_mainstems |>
filter(str_detect(habitat, "spawning")) |>
group_by(river_group) |>
summarize() |>
mutate(length_ft = st_length(geometry) |> units::set_units("ft") |> units::drop_units())
cv_mainstems_rearing <-
habistat::cv_mainstems |>
filter(str_detect(habitat, "rearing")) |>
group_by(river_group) |>
summarize() |>
mutate(length_ft = st_length(geometry) |> units::set_units("ft") |> units::drop_units())
cvpia_reach_lengths <-
bind_rows("spawn" = st_drop_geometry(cv_mainstems_spawning),
"juv" = st_drop_geometry(cv_mainstems_rearing),
"floodplain" = st_drop_geometry(cv_mainstems_rearing),
.id = "hab") |>
mutate(hab = factor(hab, levels = c("spawn", "juv", "floodplain")))
dsm_habitat_combined <-
bind_rows("instream_wua" =
dsm_habitat_instream |>
filter(hab %in% c("spawn", "juv")) |>
drop_na(wua_per_lf),
"instream_sqm" =
dsm_habitat_instream_sqm |>
filter(hab %in% c("spawn", "juv")) |>
drop_na(suitable_ac),
"floodplain_acres" =
dsm_habitat_floodplain |>
select(river_group, hab, run, flow_cfs, suitable_ac) |>
drop_na(suitable_ac),
.id = "source_datatype") |>
inner_join(cvpia_reach_lengths, by=join_by(river_group, hab)) |>
arrange(river_group, hab, run, flow_cfs) |>
group_by(river_group, hab, run) |>
mutate(wua_per_lf = zoo::na.approx(wua_per_lf, x = flow_cfs, na.rm=F),
suitable_ac = zoo::na.approx(suitable_ac, x = flow_cfs, na.rm=F)) |>
mutate(suitable_ac = coalesce(suitable_ac, (wua_per_lf * length_ft / 43560)),
wua_per_lf = coalesce(wua_per_lf, (suitable_ac * 43560 / length_ft))) |>
mutate(regional_approx =
((hab == "floodplain") & ((river_group %in% c(deer_creek_fp_proxy, cottonwood_creek_fp_proxy, tuolumne_river_fp_proxy)))) |
((hab == "juv") & ((river_group %in% c(regional_approx_groups, regional_approx_groups_sanjoaquin)))) |
((hab == "spawn") & ((river_group %in% c(regional_approx_groups_spawning, regional_approx_groups_sanjoaquin)))))|>
mutate(hab_type = case_when(hab == "juv" ~ "rearing",
hab == "spawn" ~ "spawning",
hab == "floodplain" ~ "rearing"),
hab_subtype = case_when(hab == "juv" ~ "instream",
hab == "spawn" ~ "instream",
hab == "floodplain" ~ "floodplain"))
wua_predicted_cv_mainstems_grouped |>
filter(habitat == "rearing") |>
filter(river_group %in% dsm_habitat_combined$river_group) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_ribbon(aes(ymin = wua_acres_pred_SD, ymax = wua_acres_pred_SN,
fill = paste("habistat", habitat)), alpha=0.33) +
geom_line(aes(y = wua_acres_pred, color = paste("habistat", habitat))) +
scale_x_log10() +
geom_line(data=dsm_habitat_combined |>
filter(run == "fall") |>
filter(hab_type == "rearing") |>
filter(flow_cfs <= 15000),
aes(x = flow_cfs, y = suitable_ac, color = paste("DSMhabitat", hab_subtype), linetype = regional_approx)) +
theme(legend.position = "top", panel.grid.minor = element_blank()) +
scale_color_brewer(name = "Habitat Data Source", palette = "Dark2", aesthetics = c("color", "fill")) +
scale_linetype_manual(name = "Regional Approx.", values = c("TRUE" = "dashed", "FALSE" = "solid")) +
ylab("Predicted Total Habitat (acres)") + xlab("Flow at Outlet (cfs)") +
guides(color = guide_legend(nrow = 2), linetype = guide_legend(nrow = 2)) +
ggtitle("Rearing Comparison - Total Acres")
wua_predicted_cv_mainstems_grouped |>
filter(habitat == "spawning") |>
filter(river_group %in% dsm_habitat_combined$river_group) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_ribbon(aes(ymin = wua_acres_pred_SD, ymax = wua_acres_pred_SN,
fill = paste("habistat", habitat)), alpha=0.33) +
geom_line(aes(y = wua_acres_pred, color = paste("habistat", habitat))) +
scale_x_log10() +
geom_line(data=dsm_habitat_combined |>
filter(run == "fall") |>
filter(hab_type == "spawning") |>
filter(flow_cfs <= 15000),
aes(x = flow_cfs, y = suitable_ac, color = paste("DSMhabitat", "spawning"), linetype = regional_approx)) +
theme(legend.position = "top", panel.grid.minor = element_blank()) +
scale_color_brewer(name = "Habitat Data Source", palette = "Dark2", aesthetics = c("color", "fill")) +
scale_linetype_manual(name = "Regional Approx.", values = c("TRUE" = "dashed", "FALSE" = "solid")) +
ylab("Predicted Total Habitat (acres)") + xlab("Flow at Outlet (cfs)") +
guides(color = guide_legend(nrow = 2), linetype = guide_legend(nrow = 2)) +
ggtitle("Spawning Comparison - Total Acres")
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
wua_predicted_cv_mainstems_grouped |>
filter(habitat == "rearing") |>
filter(river_group %in% dsm_habitat_combined$river_group) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_ribbon(aes(ymin = wua_per_lf_pred_SD, ymax = wua_per_lf_pred_SN,
fill = paste("habistat", habitat)), alpha=0.33) +
geom_line(aes(y = wua_per_lf_pred, color = paste("habistat", habitat))) +
scale_x_log10() +
geom_line(data=dsm_habitat_combined |>
filter(run == "fall") |>
filter(hab_type == "rearing") |>
filter(flow_cfs <= 15000),
aes(x = flow_cfs, y = wua_per_lf, color = paste("DSMhabitat", hab_subtype), linetype = regional_approx)) +
theme(legend.position = "top", panel.grid.minor = element_blank()) +
scale_color_brewer(name = "Habitat Data Source", palette = "Dark2", aesthetics = c("color", "fill")) +
scale_linetype_manual(name = "Regional Approx.", values = c("TRUE" = "dashed", "FALSE" = "solid")) +
ylab("Predicted Habitat (ft2) per LF") + xlab("Flow at Outlet (cfs)") +
guides(color = guide_legend(nrow = 2), linetype = guide_legend(nrow = 2)) +
ggtitle("Rearing Comparison - Suitable Area per Linear Ft")
wua_predicted_cv_mainstems_grouped |>
filter(habitat == "spawning") |>
filter(river_group %in% dsm_habitat_combined$river_group) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_ribbon(aes(ymin = wua_per_lf_pred_SD, ymax = wua_per_lf_pred_SN,
fill = paste("habistat", habitat)), alpha=0.33) +
geom_line(aes(y = wua_per_lf_pred, color = paste("habistat", habitat))) +
scale_x_log10() +
geom_line(data=dsm_habitat_combined |>
filter(run == "fall") |>
filter(hab_type == "spawning") |>
filter(flow_cfs <= 15000),
aes(x = flow_cfs, y = wua_per_lf, color = paste("DSMhabitat", "spawning"), linetype = regional_approx)) +
theme(legend.position = "top", panel.grid.minor = element_blank()) +
scale_color_brewer(name = "Habitat Data Source", palette = "Dark2", aesthetics = c("color", "fill")) +
scale_linetype_manual(name = "Regional Approx.", values = c("TRUE" = "dashed", "FALSE" = "solid")) +
ylab("Predicted Habitat (ft2) per LF") + xlab("Flow at Outlet (cfs)") +
guides(color = guide_legend(nrow = 2), linetype = guide_legend(nrow = 2)) +
ggtitle("Spawning Comparison - Suitable Area per Linear Ft")
## Warning in scale_x_log10(): log-10 transformation introduced infinite values.
upper_mid_sac_tribs <- DSMhabitat::modeling_exist %>%
left_join(DSMhabitat::watershed_regions, by = c('Watershed' = 'watershed')) %>%
filter(region == 'Upper-mid Sacramento River',
Watershed != 'Upper-mid Sacramento River'
)
watersheds_without_modeling <- upper_mid_sac_tribs %>%
select(watershed = Watershed, contains('spawn')) %>%
gather(species, spawn_modeling, -watershed) %>%
group_by(watershed) %>%
summarise(no_modeling = sum(spawn_modeling, na.rm = T) == 0) %>%
filter(no_modeling) %>%
pull(watershed)
cat(c("Watersheds without any Spawning (spawn) Model:", paste(watersheds_without_modeling, collapse = ", "), "\n"))
## Watersheds without any Spawning (spawn) Model: Antelope Creek, Bear Creek, Big Chico Creek, Cow Creek, Elder Creek, Mill Creek, Paynes Creek, Stony Creek, Thomes Creek
watersheds_with_spawn <- upper_mid_sac_tribs %>%
filter(!(Watershed %in% watersheds_without_modeling)) %>%
pull(Watershed)
cat(c("Watersheds with Spawning (FR spawn) Model:", paste(watersheds_with_spawn, collapse = ", "), "\n"))
## Watersheds with Spawning (FR spawn) Model: Battle Creek, Butte Creek, Clear Creek, Cottonwood Creek, Deer Creek
watersheds_with_fry <- upper_mid_sac_tribs %>%
filter(FR_fry) %>%
pull(Watershed)
cat(c("Watersheds with Rearing (FR fry) Model:", paste(watersheds_with_fry, collapse = ", "), "\n"))
## Watersheds with Rearing (FR fry) Model: Battle Creek, Clear Creek, Cottonwood Creek, Cow Creek
watersheds_with_juv <- upper_mid_sac_tribs %>%
filter(FR_juv) %>%
pull(Watershed)
cat(c("Watersheds with Rearing (FR juv) Model:", paste(watersheds_with_juv, collapse = ", "), "\n"))
## Watersheds with Rearing (FR juv) Model: Battle Creek, Clear Creek, Cottonwood Creek, Cow Creek, Deer Creek
river_lengths <-
mainstems_comid |> # note this is a general, not run-specific habitat extent dataset
st_drop_geometry() |>
group_by(watershed, habitat) |>
summarize(length_ft = sum(length_ft)) |>
pivot_wider(names_from = habitat,
values_from = length_ft,
names_repair = janitor::make_clean_names) |>
transmute(river_group = watershed,
juv = coalesce(rearing, 0) + coalesce(rearing_and_spawning, 0),
spawn = coalesce(rearing_and_spawning, 0)) |>
pivot_longer(cols = c(juv, spawn),
names_to = "hab",
values_to = "length_ft") |>
mutate(length_rm = length_ft / 5280) |>
glimpse()
## `summarise()` has grouped output by 'watershed'. You can override using the
## `.groups` argument.
## Rows: 48
## Columns: 5
## Groups: watershed [24]
## $ watershed <chr> "American River", "American River", "Antelope Creek", "Ant…
## $ river_group <chr> "American River", "American River", "Antelope Creek", "Ant…
## $ hab <chr> "juv", "spawn", "juv", "spawn", "juv", "spawn", "juv", "sp…
## $ length_ft <dbl> 118816.94, 99546.05, 160133.24, 126501.58, 127772.02, 1215…
## $ length_rm <dbl> 22.50321, 18.85342, 30.32826, 23.95863, 24.19925, 23.02454…
interp_flows <- seq_log10(50, 15000, 0.05, snap=100) # to align with habistat flows
regional_approx_curve <-
DSMhabitat::upper_mid_sac_region_instream |>
# select spawn and juv and convert WUA/1000LF to WUA/LF
transmute(flow_cfs,
spawn = FR_spawn_wua / 1000,
juv = FR_juv_wua / 1000) |>
pivot_longer(cols = c(spawn, juv), names_to = "hab", values_to = "wua_per_lf") |>
# interpolate onto the habistat flow sequence
group_by(hab) |>
reframe(wua_per_lf = approx(y = wua_per_lf, x = log(flow_cfs), xout = log(interp_flows), na.rm = F)$y,
flow_cfs = interp_flows) |>
drop_na()
regional_approx_curve |>
ggplot() +
geom_line(aes(x = flow_cfs, y = wua_per_lf, color = hab)) +
scale_x_log10() +
scale_y_continuous(limits = c(0, NA),
breaks = scales::breaks_width(2),
expand = c(0, 0)) +
annotation_logticks(sides = "b") +
theme(panel.grid.minor = element_blank(),
strip.placement = "outside",
strip.text.y.left = element_text(angle = 0)) +
labs(title = "Upper-Mid Sacramento Regional Approximation",
subtitle = "Fall Run Chinook Salmon",
caption = paste("Spawning (spawn) based on: Battle Creek, Butte Creek, Clear Creek.",
"Rearing (juv) based on: Battle Creek, Butte Creek, Clear Creek, Cow Creek.",
sep = "\n")) +
xlab("Flow (cfs)") +
ylab("Suitable Habitat Area (ft2) per linear ft") +
scale_color_brewer(name = "Habitat Type", palette = "Dark2", aesthetics = c("color", "fill"))
instream_regional_approx_est <-
river_lengths |>
inner_join(regional_approx_curve,
by = join_by(hab),
relationship = "many-to-many") |>
mutate(suitable_ac = wua_per_lf * length_ft / 43560,
habitat = case_when(hab == "spawn" ~ "spawning",
hab == "juv" ~ "rearing"))
instream_regional_approx_est |>
filter(river_group %in% dsm_habitat_combined$river_group) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_line(aes(y = suitable_ac, color = habitat)) +
scale_x_log10() +
theme(legend.position = "top", panel.grid.minor = element_blank()) +
scale_color_brewer(name = "Habitat Type", palette = "Dark2", aesthetics = c("color", "fill")) +
scale_linetype_manual(name = "Regional Approx.", values = c("TRUE" = "dashed", "FALSE" = "solid")) +
ylab("Predicted Total Habitat (acres)") + xlab("Flow at Outlet (cfs)") +
guides(color = guide_legend(nrow = 2), linetype = guide_legend(nrow = 2)) +
ggtitle("Instream Regional Approximation Applied to All Watersheds")
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's linetype values.
wua_predicted_cv_mainstems_grouped |>
filter(habitat == "rearing") |>
filter(river_group %in% dsm_habitat_combined$river_group) |>
filter(river_group %in% upper_mid_sac_tribs$Watershed) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_ribbon(aes(ymin = wua_acres_pred_SD, ymax = wua_acres_pred_SN,
fill = "habistat"), alpha=0.33) +
geom_line(aes(y = wua_acres_pred,
linetype = "combined instream + floodplain",
color = "habistat")) +
geom_line(data=dsm_habitat_combined |>
filter(run == "fall") |>
filter(hab_type == "rearing") |>
filter(flow_cfs <= 15000) |>
#filter(!regional_approx) |>
filter(river_group %in% upper_mid_sac_tribs$Watershed),
aes(x = flow_cfs,
y = suitable_ac,
linetype = hab_subtype,
color = paste("DSMHabitat", if_else(regional_approx, "approximation", "hydraulic model")))) +
geom_line(data = instream_regional_approx_est |>
filter(habitat == "rearing") |>
filter(river_group %in% upper_mid_sac_tribs$Watershed),
aes(y = suitable_ac,,
linetype = "instream",
color = "DSMHabitat regional approx")) +
scale_x_log10() +
scale_linetype_manual(name = "Habitat Type",
values = c("instream" = "solid", #"#6388b4",
"floodplain" = "dashed", #"#ef6f6a",
"combined instream + floodplain" = "dotted")) + ##bb7693"),) +
scale_color_manual(name = "Data Source",
values = c("habistat" = "#ffae34",
"DSMHabitat hydraulic model" = "#6388b4",
"DSMHabitat regional approx" = "#ef6f6a"),
aesthetics = c("color", "fill")) +
theme(legend.position = "top",
panel.grid.minor = element_blank()) +
guides(color = guide_legend(nrow = 3),
linetype = guide_legend(nrow = 3)) +
ylab("Predicted Total Habitat (acres)") +
xlab("Flow at Outlet (cfs)") +
ggtitle("Rearing Comparison for Upper-Mid Sacramento Tributaries - Total Acres")
wua_predicted_cv_mainstems_grouped |>
filter(habitat == "spawning") |>
filter(river_group %in% dsm_habitat_combined$river_group) |>
filter(river_group %in% upper_mid_sac_tribs$Watershed) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_ribbon(aes(ymin = wua_acres_pred_SD, ymax = wua_acres_pred_SN,
fill = "habistat"), alpha=0.33) +
geom_line(aes(y = wua_acres_pred,
linetype = "instream",
color = "habistat")) +
geom_line(data=dsm_habitat_combined |>
filter(run == "fall") |>
filter(hab_type == "spawning") |>
filter(flow_cfs <= 15000) |>
#filter(!regional_approx) |>
filter(river_group %in% upper_mid_sac_tribs$Watershed),
aes(x = flow_cfs,
y = suitable_ac,
linetype = hab_subtype,
color = paste("DSMHabitat", if_else(regional_approx, "approximation", "hydraulic model")))) +
geom_line(data = instream_regional_approx_est |>
filter(habitat == "spawning") |>
filter(river_group %in% upper_mid_sac_tribs$Watershed),
aes(y = suitable_ac,,
linetype = "instream",
color = "DSMHabitat regional approx")) +
scale_x_log10() +
scale_linetype_manual(name = "Habitat Type",
values = c("instream" = "solid", #"#6388b4",
"floodplain" = "dashed", #"#ef6f6a",
"combined instream + floodplain" = "dotted")) + ##bb7693"),) +
scale_color_manual(name = "Data Source",
values = c("habistat" = "#ffae34",
"DSMHabitat hydraulic model" = "#6388b4",
"DSMHabitat regional approx" = "#ef6f6a"),
aesthetics = c("color", "fill")) +
theme(legend.position = "top",
panel.grid.minor = element_blank()) +
guides(color = guide_legend(nrow = 3),
linetype = guide_legend(nrow = 3)) +
ylab("Predicted Total Habitat (acres)") +
xlab("Flow at Outlet (cfs)") +
ggtitle("Spawning Comparison for Upper-Mid Sacramento Tributaries - Total Acres")
# from DSMhabitat data-raw/watershed/floodplain_utils.R and modified for this use case
# attempt at dec_jun_mean_flow_scaling value replication...
# DSMflow::flows_cfs$biop_itp_2018_2019 |> select(date, `Antelope Creek`, `Deer Creek`) |> janitor::clean_names() |> filter(month(date) %in% c(6,12)) |> mutate(scale = antelope_creek/deer_creek)
scale_fp_flow_area <- function(ws, proxy_ws, flow_scalar) {
watershed_metadata <- filter(DSMhabitat::floodplain_modeling_metadata,
watershed == ws)
# optional manual overrides of values that are pulled from watershed_metadata:
selected_proxy_ws <- if(missing(proxy_ws)) watershed_metadata$scaling_watershed else proxy_ws
selected_flow_scalar <- if(missing(flow_scalar)) watershed_metadata$dec_jun_mean_flow_scaling else flow_scalar
low_gradient <- filter(DSMhabitat::low_gradient_lengths, watershed_name == ws)
rearing_extents <- filter(DSMhabitat::watershed_lengths, watershed == ws, lifestage == "rearing") %>%
select(watershed, species, miles) %>%
spread(species, miles)
proxy_watershed_metadata <- filter(DSMhabitat::floodplain_modeling_metadata,
watershed == selected_proxy_ws)
species_present <- filter(DSMhabitat::watershed_species_present,
watershed_name == ws)
proxy_data <- switch(selected_proxy_ws,
"Deer Creek" = DSMhabitat::deer_creek_floodplain,
"Tuolumne River" = DSMhabitat::tuolumne_river_floodplain,
"Cottonwood Creek" = DSMhabitat::cottonwood_creek_floodplain)
# scale flow
scaled_flow <- round(proxy_data$flow_cfs * selected_flow_scalar)
# fall run area
# divide floodplain area by watershed length of proxy watershed to get area/mile, scale to hydrology
scaled_area_per_mile_FR <- (proxy_data$FR_floodplain_acres / proxy_watershed_metadata$modeled_length_mi) *
selected_flow_scalar
# apportion area by high gradient/low gradient, .1 is downscaling for high gradient
fp_area_FR <- round((scaled_area_per_mile_FR * low_gradient$fr) +
(scaled_area_per_mile_FR * (rearing_extents$fr - low_gradient$fr) * 0.1), 2)
result <- data.frame(
flow_cfs = scaled_flow,
FR_floodplain_acres = fp_area_FR
)
if (species_present$lfr) {
# latefall floodplain area
lfr_proxy <- if ("LFR_floodplain_acres" %in% colnames(proxy_data)) {
proxy_data$LFR_floodplain_acres
} else {
proxy_data$FR_floodplain_acres
}
scaled_area_per_mile_LFR <- (lfr_proxy / proxy_watershed_metadata$modeled_length_mi) *
selected_flow_scalar
fp_area_LFR <-round((scaled_area_per_mile_LFR * low_gradient$lfr) +
(scaled_area_per_mile_LFR * (rearing_extents$lfr - low_gradient$lfr) * 0.1), 2)
result <- bind_cols(result, LFR_floodplain_acres = fp_area_LFR)
}
if (species_present$wr) {
# winter run floodplain area
wr_proxy <- if ("WR_floodplain_acres" %in% colnames(proxy_data)) {
proxy_data$WR_floodplain_acres
} else {
proxy_data$FR_floodplain_acres
}
scaled_area_per_mile_WR <- (wr_proxy / proxy_watershed_metadata$modeled_length_mi) *
selected_flow_scalar
fp_area_WR <-round((scaled_area_per_mile_WR * low_gradient$wr) +
(scaled_area_per_mile_WR * (rearing_extents$wr - low_gradient$wr) * 0.1), 2)
result <- bind_cols(result, WR_floodplain_acres = fp_area_WR)
}
if (species_present$sr) {
# spring run floodplain area
sr_proxy <- if ("SR_floodplain_acres" %in% colnames(proxy_data)) {
proxy_data$SR_floodplain_acres
} else {
if ("ST_floodplain_acres" %in% colnames(proxy_data)) {
proxy_data$ST_floodplain_acres
} else {
proxy_data$FR_floodplain_acres
}
}
scaled_area_per_mile_SR <- (sr_proxy / proxy_watershed_metadata$modeled_length_mi) *
selected_flow_scalar
fp_area_SR <- round((scaled_area_per_mile_SR * low_gradient$sr) +
(scaled_area_per_mile_SR * (rearing_extents$sr - low_gradient$sr) * 0.1), 2)
result <- bind_cols(result, SR_floodplain_acres = fp_area_SR)
}
if (species_present$st) {
# steelhead floodplain area
st_proxy <- if ("ST_floodplain_acres" %in% colnames(proxy_data)) {
proxy_data$ST_floodplain_acres
} else {
if ("SR_floodplain_acres" %in% colnames(proxy_data)) {
proxy_data$SR_floodplain_acres
} else {
proxy_data$FR_floodplain_acres
}
}
scaled_area_per_mile_ST <- (st_proxy / proxy_watershed_metadata$modeled_length_mi) *
selected_flow_scalar
fp_area_ST <- round((scaled_area_per_mile_ST * low_gradient$st) +
(scaled_area_per_mile_ST * (rearing_extents$st - low_gradient$st) * 0.1), 2)
result <- bind_cols(result, ST_floodplain_acres = fp_area_ST)
}
result <- result |>
mutate(watershed = ws,
proxy_watershed = selected_proxy_ws)
return(mutate(result, watershed = ws))
}
# comparison file:
# wua_predicted_cv_mainstems_grouped
watersheds_fp_scale <- DSMhabitat::floodplain_modeling_metadata |>
filter(!is.na(dec_jun_mean_flow_scaling)) |>
pull(watershed)
full_fp_results <- data.frame()
for(i in 1:length(watersheds_fp_scale)) {
result <- scale_fp_flow_area(watersheds_fp_scale[i])
full_fp_results <- bind_rows(full_fp_results, result)
}
full_fp_results |>
ggplot(aes(x = flow_cfs, y = FR_floodplain_acres, color = proxy_watershed)) +
geom_line() +
facet_wrap(~watershed, scales = "free_y") +
scale_x_log10(breaks = scales::breaks_log(),
labels = scales::label_number()) +
theme(panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 90))
# full_deer_creek_fp_results <- data.frame()
# for(i in 1:length(deer_creek_fp_proxy)) {
# result <- scale_fp_flow_area(deer_creek_fp_proxy[i])
# full_deer_creek_fp_results <- bind_rows(full_deer_creek_fp_results, result)
# }
#
# flows <- c(DSMhabitat::mill_creek_floodplain$flow_cfs, DSMhabitat::cow_creek_floodplain$flow_cfs)
#
# # get_approx <- function(df, species_column) {
# # approxfun(df$flow_cfs, df[,species_column, drop = TRUE], rule = 2)
# # }
#
# full_deer_creek_fp_results |>
# select(flow_cfs, floodplain_acres = FR_floodplain_acres, watershed) |>
# ggplot(aes(x = flow_cfs, y = floodplain_acres, color = watershed)) +
# geom_line() +
# ggtitle("Deer Creek Floodplain Proxy Approximations")
#
# full_cottonwood_creek_fp_results <- data.frame()
# for(i in 1:length(cottonwood_creek_fp_proxy)) {
# result <- scale_fp_flow_area(cottonwood_creek_fp_proxy[i])
# full_cottonwood_creek_fp_results <- bind_rows(full_cottonwood_creek_fp_results, result)
# }
#
# full_cottonwood_creek_fp_results |>
# select(flow_cfs, floodplain_acres = FR_floodplain_acres, watershed) |>
# ggplot(aes(x = flow_cfs, y = floodplain_acres, color = watershed)) +
# geom_line() +
# ggtitle("Cottonwood Creek Floodplain Proxy Approximations")
floodplain_proxy_est_v1 <-
DSMhabitat::floodplain_modeling_metadata |>
filter(!is.na(dec_jun_mean_flow_scaling)) |>
mutate(result = pmap(list(watershed, scaling_watershed, dec_jun_mean_flow_scaling),
scale_fp_flow_area)) |>
rename(river_group = watershed) |>
unnest(result) |>
mutate(habitat = "floodplain") |>
pivot_longer(ends_with("_floodplain_acres"),
names_transform = \(x) str_replace(x, "_floodplain_acres", ""),
names_to = "run",
values_to = "floodplain_ac") |>
mutate(run = factor(run,
levels = c("FR", "LFR", "WR", "SR", "ST"),
labels = c("fall", "late fall", "winter", "spring", "steelhead"))) |>
mutate(suitable_ac = #if_else(river_group %in% suitability_already_applied,
floodplain_ac,
# DSMhabitat::apply_suitability(floodplain_ac * 4046.86) / 4046.86)
)
wua_predicted_cv_mainstems_grouped |>
filter(habitat == "rearing") |>
filter(river_group %in% dsm_habitat_combined$river_group) |>
filter(river_group %in% upper_mid_sac_tribs$Watershed) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_ribbon(aes(ymin = wua_acres_pred_SD, ymax = wua_acres_pred_SN,
fill = "habistat"), alpha=0.33) +
geom_line(aes(y = wua_acres_pred,
linetype = "combined instream + floodplain",
color = "habistat")) +
geom_line(data=dsm_habitat_combined |>
filter(run == "fall") |>
filter(hab_type == "rearing") |>
filter(flow_cfs <= 15000) |>
#filter(!regional_approx) |>
filter(river_group %in% upper_mid_sac_tribs$Watershed),
aes(x = flow_cfs,
y = suitable_ac,
linetype = hab_subtype,
color = paste("DSMHabitat", if_else(regional_approx, "approximation", "hydraulic model")))) +
geom_line(data = instream_regional_approx_est |>
filter(habitat == "rearing") |>
filter(river_group %in% upper_mid_sac_tribs$Watershed),
aes(y = suitable_ac,,
linetype = "instream",
color = "DSMHabitat approximation")) +
geom_line(data = floodplain_proxy_est_v1 |>
filter(habitat == "floodplain" & run == "fall") |>
filter(river_group %in% upper_mid_sac_tribs$Watershed),
aes(y = suitable_ac,,
linetype = "floodplain",
color = "DSMHabitat approximation")) +
scale_x_log10() +
scale_linetype_manual(name = "Habitat Type",
values = c("instream" = "solid", #"#6388b4",
"floodplain" = "dashed", #"#ef6f6a",
"combined instream + floodplain" = "dotted")) + ##bb7693"),) +
scale_color_manual(name = "Data Source",
values = c("habistat" = "#ffae34",
"DSMHabitat hydraulic model" = "#6388b4",
"DSMHabitat approximation" = "#ef6f6a"),
aesthetics = c("color", "fill")) +
theme(legend.position = "top",
panel.grid.minor = element_blank()) +
guides(color = guide_legend(nrow = 3),
linetype = guide_legend(nrow = 3)) +
ylab("Predicted Total Habitat (acres)") +
xlab("Flow at Outlet (cfs)") +
ggtitle("Rearing Comparison for Upper-Mid Sacramento Tributaries - Total Acres")
Replicate monthly flow scaling
Streams: https://www.sciencebase.gov/catalog/item/61e7442ed34e3618e01cf68f Monthly inflows: https://www.sciencebase.gov/catalog/item/61e74499d34e3618e01cf694 Monthly diversions: https://www.sciencebase.gov/catalog/item/61e74461d34e3618e01cf691
Version using CVHM2:
inflow_locs <-
read_sf(here::here("data-raw", "source", "monthly_flows", "CVHM2_InflowLocations.shp.zip")) |>
janitor::clean_names()
inflow_vals <-
read_csv(here::here("data-raw", "source", "monthly_flows", "CVHM2_Inflows_1921_2019.csv")) |>
pivot_longer(cols = c(-Month, -Year, -YYYYMM),
names_to = "inflow_id",
values_to = "flow_cmd") |>
janitor::clean_names()
## Rows: 1181 Columns: 65
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (65): Month, Year, YYYYMM, SACR_205, COWC_211, BATT_220, COTT_NF, COTT_M...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
inflow_xw <-
read_csv(here::here("data-raw", "source", "monthly_flows", "CVHM2_Inflows_XW.csv"))
## Rows: 65 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): name, inflow_id, river_group
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mean_monthly_inflow <-
inflow_vals |>
group_by(inflow_id, month) |>
summarize(flow_cmd = mean(flow_cmd)) |>
# convert cubic meters per day to cubic feet per second
mutate(flow_cfs = flow_cmd * (1 / 0.3048)^3 * (1 / 86400)) |>
inner_join(inflow_xw, by = join_by(inflow_id)) %>%
inner_join(. |>
filter(inflow_id %in% c("DEER_256", "COTT_MF", "TUOL_135")) |>
select(inflow_id, month, flow_cfs) |>
pivot_wider(names_from = inflow_id, values_from = flow_cfs) |>
rename(q_deer = DEER_256,
q_cottonwood = COTT_MF,
q_tuolumne = TUOL_135),
by = join_by(month))
## `summarise()` has grouped output by 'inflow_id'. You can override using the
## `.groups` argument.
flow_scalars <-
mean_monthly_inflow |>
filter(month %in% c(12, 1, 2, 3, 4, 5, 6)) |>
group_by(river_group) |>
summarize(across(c(flow_cfs, starts_with("q_")), mean)) |>
mutate(qscal_cottonwood = flow_cfs / q_cottonwood,
qscal_deer = flow_cfs / q_deer,
qscal_tuolumne = flow_cfs / q_tuolumne)
# flow_scalars |>
# rename(scaling_flow = flow_cfs) |>
# pivot_longer(cols = starts_with("qscal_"),
# names_transform = \(x) str_replace(x, "qscal_", ""),
# values_to = "flow_scalar") |>
# mutate(proxy_ws = case_when(name == "cottonwood" ~ "Cottonwood Creek",
# name == "deer" ~ "Deer Creek",
# name == "tuolumne" ~ "Tuolumne River")) |>
# mutate(result = pmap(list(river_group, proxy_ws, flow_scalar),
# possibly(scale_fp_flow_area, NA))) |>
# unnest(result) |>
# mutate(habitat = "floodplain") |>
# pivot_longer(ends_with("_floodplain_acres"),
# names_transform = \(x) str_replace(x, "_floodplain_acres", ""),
# names_to = "run",
# values_to = "suitable_ac") |>
# mutate(run = factor(run,
# levels = c("FR", "LFR", "WR", "SR", "ST"),
# labels = c("fall", "late fall", "winter", "spring", "steelhead")))
Version using DSMFlow:
mean_dec_jun_flows <-
DSMflow::flows_cfs$biop_itp_2018_2019 |>
pivot_longer(cols = -date, names_to = "watershed", values_to = "flow_cfs") |>
mutate(water_year = year(date) + if_else(month(date) >= 10, 1, 0)) |>
filter(month(date) %in% c(12, 1, 2, 3, 4, 5, 6)) |>
group_by(watershed) |>
summarize(mean_flow_cfs = mean(flow_cfs))
# select(date, `Antelope Creek`, `Deer Creek`) |>
# janitor::clean_names() |>
# filter(month(date) %in% c(6,12)) |> mutate(scale = antelope_creek/deer_creek)
proxy_watersheds <-
c("dc" = "Deer Creek",
"cc" = "Cottonwood Creek",
"tr" = "Tuolumne River")
proxy_flows <-
mean_dec_jun_flows |>
filter(watershed %in% proxy_watersheds) |>
pivot_wider(names_from = watershed,
values_from = mean_flow_cfs) |>
rename_with(\(x) setNames(names(proxy_watersheds), proxy_watersheds)[x])
flow_scalars <-
mean_dec_jun_flows |>
mutate(qscal_dc = mean_flow_cfs / proxy_flows$dc,
qscal_cc = mean_flow_cfs / proxy_flows$cc,
qscal_tr = mean_flow_cfs / proxy_flows$tr) |>
left_join(DSMhabitat::watershed_regions, by = join_by(watershed)) |>
left_join(DSMhabitat::floodplain_modeling_metadata |>
select(watershed, scaling_watershed, dec_jun_mean_flow_scaling),
by = join_by(watershed)) |>
mutate(selected_proxy = case_when(
watershed %in% proxy_watersheds ~ watershed,
!is.na(scaling_watershed) ~ scaling_watershed,
region %in% c("South Delta", "San Joaquin River") ~ "Tuolumne River",
abs(mean_flow_cfs - proxy_flows$dc) <= abs(mean_flow_cfs - proxy_flows$cc) ~ "Deer Creek",
abs(mean_flow_cfs - proxy_flows$dc) > abs(mean_flow_cfs - proxy_flows$cc) ~ "Cottonwood Creek"))
proxy_lookup <-
flow_scalars |>
select(watershed, selected_proxy) |>
deframe()
floodplain_proxy_est <-
flow_scalars |>
rename(river_group = watershed) |>
filter(!str_detect(river_group, "Sacramento") & !str_detect(river_group, "San Joaquin")) |>
pivot_longer(cols = starts_with("qscal_"),
names_transform = \(x) proxy_watersheds[str_replace(x, "qscal_", "")],
names_to = "proxy_ws",
values_to = "flow_scalar") |>
mutate(selected = (proxy_ws == proxy_lookup[river_group]),
flow_scalar_combined = if_else(selected,
coalesce(dec_jun_mean_flow_scaling, flow_scalar),
flow_scalar),
result = pmap(list(river_group, proxy_ws, flow_scalar_combined),
possibly(scale_fp_flow_area, NA))) |>
drop_na(result) |>
unnest(result)
floodplain_proxy_est |>
ggplot() +
facet_wrap(~watershed_labels[river_group], scales = "free_y") +
geom_line(aes(x = flow_cfs, y = FR_floodplain_acres, color = proxy_ws)) +
scale_x_log10(breaks = scales::breaks_log(),
labels = scales::label_number()) +
theme(panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 90))
floodplain_proxy_est |>
filter(selected) |>
ggplot() +
facet_wrap(~watershed_labels[river_group], scales = "free_y") +
geom_line(aes(x = flow_cfs, y = FR_floodplain_acres, color = proxy_ws)) +
scale_x_log10(breaks = scales::breaks_log(),
labels = scales::label_number()) +
theme(panel.grid.minor = element_blank(),
axis.text.x = element_text(angle = 90))
floodplain_proxy_est_pivot <-
floodplain_proxy_est |>
filter(selected) |>
mutate(habitat = "floodplain") |>
pivot_longer(ends_with("_floodplain_acres"),
names_transform = \(x) str_replace(x, "_floodplain_acres", ""),
names_to = "run",
values_to = "floodplain_ac") |>
mutate(run = factor(run,
levels = c("FR", "LFR", "WR", "SR", "ST"),
labels = c("fall", "late fall", "winter", "spring", "steelhead"))) |>
mutate(suitable_ac = if_else(river_group %in% suitability_already_applied,
floodplain_ac,
DSMhabitat::apply_suitability(floodplain_ac * 4046.86) / 4046.86))
# version summing instream and floodplain
dsm_habitat_combined_sum <-
dsm_habitat_combined |>
ungroup() |>
filter(hab_type == "rearing" & hab %in% c("juv", "floodplain") & !regional_approx) |> # only produces this combined estimate when there are modeled values for both instream AND floodplain
select(river_group, run, hab, flow_cfs, wua_per_lf, suitable_ac) |>
group_by(river_group, run) |>
arrange(river_group, run, hab, flow_cfs) |>
complete(hab = c("juv", "floodplain"), flow_cfs) |>
complete(hab = c("juv", "floodplain"), flow_cfs = interp_flows) |>
group_by(river_group, run, hab) |>
mutate(across(c(wua_per_lf, suitable_ac),
\(v) zoo::na.approx(v, x = log(flow_cfs), na.rm = F, rule = 2:2))) |>
group_by(river_group, run, flow_cfs) |>
summarize(across(c(wua_per_lf, suitable_ac), sum))
## `summarise()` has grouped output by 'river_group', 'run'. You can override
## using the `.groups` argument.
dsm_habitat_combined_sum |>
filter(run == "fall") |>
ggplot() +
facet_wrap(~watershed_labels[river_group], scales = "free_y") +
geom_line(data = dsm_habitat_combined |>
filter(run == "fall" & hab %in% c("juv", "floodplain")),
aes(x = flow_cfs, y = suitable_ac, linetype = "original", color = hab)) +
geom_line(aes(x = flow_cfs, y = suitable_ac, linetype = "combined")) +
scale_x_log10() +
scale_linetype_manual(values = c("original" = "solid", "combined" = "dashed"))
dsm_habitat_proxy_estimates_combined <-
bind_rows("juv" = instream_regional_approx_est |>
filter(hab == "juv") |>
ungroup() |>
select(river_group, flow_cfs, suitable_ac),
"floodplain" = floodplain_proxy_est_pivot |>
ungroup() |>
filter(run == "fall") |>
select(river_group, flow_cfs, suitable_ac),
.id = "hab")
dsm_habitat_proxy_estimates_combined_sum <-
dsm_habitat_proxy_estimates_combined |>
select(river_group, hab, flow_cfs, suitable_ac) |>
group_by(river_group) |>
arrange(river_group, hab, flow_cfs) |>
complete(hab = c("juv", "floodplain"), flow_cfs) |>
complete(hab = c("juv", "floodplain"), flow_cfs = interp_flows) |>
group_by(river_group, hab) |>
mutate(suitable_ac = zoo::na.approx(suitable_ac,
x = log(flow_cfs),
na.rm = F,
rule = 2:2)) |>
group_by(river_group, flow_cfs) |>
summarize(suitable_ac = sum(suitable_ac))
## `summarise()` has grouped output by 'river_group'. You can override using the
## `.groups` argument.
dsm_habitat_proxy_estimates_combined_sum |>
ggplot() +
facet_wrap(~watershed_labels[river_group], scales = "free_y") +
geom_line(data = dsm_habitat_proxy_estimates_combined,
aes(x = flow_cfs, y = suitable_ac, linetype = "original", color = hab)) +
geom_line(aes(x = flow_cfs, y = suitable_ac, linetype = "combined")) +
scale_x_log10() +
scale_linetype_manual(values = c("original" = "solid", "combined" = "dashed"))
wua_predicted_cv_mainstems_grouped |>
filter(habitat == "rearing") |>
filter(river_group %in% dsm_habitat_combined$river_group) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free") +
geom_ribbon(aes(ymin = wua_acres_pred_SD, ymax = wua_acres_pred_SN,
fill = "habistat"), alpha=0.33) +
geom_line(aes(y = wua_acres_pred,
linetype = "combined instream + floodplain",
color = "habistat")) +
geom_line(data = instream_regional_approx_est |>
filter(habitat == "rearing"), # |>
aes(y = suitable_ac,,
linetype = "instream",
color = "DSMHabitat approximation")) +
geom_line(data = floodplain_proxy_est_pivot |>
filter(habitat == "floodplain" & run == "fall"), # |>
aes(y = suitable_ac,,
linetype = "floodplain",
color = "DSMHabitat approximation")) +
geom_line(data=dsm_habitat_combined |>
filter(run == "fall") |>
filter(hab_type == "rearing") |>
filter(flow_cfs <= 15000),
aes(x = flow_cfs,
y = suitable_ac,
linetype = hab_subtype,
color = paste("DSMHabitat", if_else(regional_approx, "approximation", "hydraulic model")))) +
geom_line(data = dsm_habitat_proxy_estimates_combined_sum,
aes(y = suitable_ac,
linetype = "combined instream + floodplain",
color = "DSMHabitat approximation")) +
scale_x_log10() +
scale_linetype_manual(name = "Habitat Type",
values = c("instream" = "solid",
"floodplain" = "dashed",
"combined instream + floodplain" = "dotted")) +
scale_color_manual(name = "Data Source",
values = c("habistat" = "#ffae34",
"DSMHabitat hydraulic model" = "#6388b4",
"DSMHabitat approximation" = "#ef6f6a"),
aesthetics = c("color", "fill")) +
theme(legend.position = "top",
panel.grid.minor = element_blank()) +
guides(color = guide_legend(nrow = 3),
linetype = guide_legend(nrow = 3)) +
ylab("Predicted Total Habitat (acres)") +
xlab("Flow at Outlet (cfs)") +
ggtitle("Rearing Comparison - Total Acres")
reach_groups <- list(
"Watersheds with DSMHabitat Instream Regional Approximation" = regional_approx_groups,
"Watersheds with DSMHabitat Floodplain Proxy Approximation" = fp_proxy_groups,
"Watersheds in Habistat training dataset" = c("Deer Creek", "Mokelumne River", "Stanislaus River", "Yuba River"),
"Watersheds with DSMHabitat full model data" = c("American River", "Bear River", "Clear Creek", "Cottonwood Creek","Feather River", "Tuolumne River"))
for (x in names(reach_groups)) {
plt <- wua_predicted_cv_mainstems_grouped |>
filter(habitat == "rearing") |>
filter(river_group %in% reach_groups[[x]]) |>
ggplot(aes(x = flow_cfs)) +
facet_wrap(~watershed_labels[river_group], scales="free_y") +
geom_line(aes(y = wua_acres_pred, color = "habistat")) +
geom_line(data = dsm_habitat_combined_sum |>
filter(run == "fall") |>
filter(river_group %in% reach_groups[[x]]) |>
filter(flow_cfs >= 50 & flow_cfs <= 15000),
aes(y = suitable_ac,
color = "DSMHabitat hydraulic model")) +
geom_line(data = dsm_habitat_proxy_estimates_combined_sum |>
filter(river_group %in% reach_groups[[x]]) |>
filter(flow_cfs >= 50 & flow_cfs <= 15000),
aes(y = suitable_ac,
color = "DSMHabitat approximation")) +
scale_color_manual(name = "Data Source",
values = c("habistat" = "#ffae34",
"DSMHabitat hydraulic model" = "#6388b4",
"DSMHabitat approximation" = "#ef6f6a"),
aesthetics = c("color", "fill")) +
scale_x_log10() +
theme(legend.position = "top",
panel.grid.minor = element_blank()) +
guides(color = guide_legend(nrow = 1)) +
ylab("Predicted Total Habitat (acres)") +
xlab("Flow at Outlet (cfs)") +
ggtitle(x)
print(plt)
}
## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 75 rows containing missing values or values outside the scale range
## (`geom_line()`).
# watersheds used in the proxy models:
ws_proxies <-
list("spawn" = c("Battle Creek", "Butte Creek", "Clear Creek",
"Calaveras River", "Mokelumne River"), # Cottonwood not used
"juv" = c("Battle Creek", "Butte Creek", "Clear Creek", "Cow Creek",
"Calaveras River", "Mokelumne River"), # Cottonwood not used
"floodplain" = c("Deer Creek", "Cottonwood Creek", "Tuolumne River"))
# watersheds that use the proxy models:
ws_uses_proxies <-
list("spawn" = c(regional_approx_groups_spawning, regional_approx_groups_sanjoaquin),
"juv" = c(regional_approx_groups, regional_approx_groups_sanjoaquin),
"floodplain" = fp_proxy_groups)
# Consumnes River derived from average WUA of Calaveras River and Mokelumne River
# watersheds that have empirical data:
ws_has_model <-
dsm_habitat_combined |>
filter(!regional_approx) |>
group_by(hab, river_group) |>
summarize(.groups = "drop") |>
nest(.by = hab) |>
mutate(data = map(data, deframe)) |>
deframe()
ws_proxy_summary <-
bind_rows("spawn" = cv_mainstems_spawning,
"juv" = cv_mainstems_rearing,
"floodplain" = cv_mainstems_rearing,
.id = "hab") |>
mutate(hab = factor(hab,
levels = c("spawn", "juv", "floodplain"),
labels = c("Spawning", "Rearing", "Floodplain"))) |>
mutate(has_model = map2_lgl(river_group, hab, \(x, y) x %in% ws_has_model[[y]]),
is_proxy = map2_lgl(river_group, hab, \(x, y) x %in% ws_proxies[[y]]),
uses_proxy = map2_lgl(river_group, hab, \(x, y) x %in% ws_uses_proxies[[y]])) |>
mutate(label = case_when(is_proxy ~ "Has Model (used as a proxy)",
has_model ~ "Has Model",
uses_proxy ~ "Does Not Have Model (uses a proxy)")) |>
drop_na(label)
ws_proxy_summary |>
st_drop_geometry() |>
select(river_group, hab, label) |>
pivot_wider(names_from = hab, values_from = label) |>
knitr::kable()
| river_group | Spawning | Rearing | Floodplain |
|---|---|---|---|
| American River | Has Model | Has Model | Has Model |
| Antelope Creek | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) |
| Battle Creek | Has Model (used as a proxy) | Has Model (used as a proxy) | Has Model |
| Bear Creek | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) |
| Bear River | Has Model | Has Model | Has Model |
| Big Chico Creek | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) | Has Model |
| Butte Creek | Has Model (used as a proxy) | Has Model (used as a proxy) | Has Model |
| Calaveras River | Has Model (used as a proxy) | Has Model (used as a proxy) | Does Not Have Model (uses a proxy) |
| Clear Creek | Has Model (used as a proxy) | Has Model (used as a proxy) | Has Model |
| Cosumnes River | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) | Has Model |
| Cottonwood Creek | Has Model | Has Model | Has Model (used as a proxy) |
| Cow Creek | Does Not Have Model (uses a proxy) | Has Model (used as a proxy) | Does Not Have Model (uses a proxy) |
| Deer Creek | Has Model | Has Model | Has Model (used as a proxy) |
| Elder Creek | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) | Has Model |
| Feather River | Has Model | Has Model | Has Model |
| Merced River | Has Model | Has Model | Has Model |
| Mill Creek | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) |
| Mokelumne River | Has Model (used as a proxy) | Has Model (used as a proxy) | Has Model |
| Paynes Creek | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) |
| Stanislaus River | Has Model | Has Model | Has Model |
| Stony Creek | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) |
| Thomes Creek | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) | Does Not Have Model (uses a proxy) |
| Tuolumne River | Has Model | Has Model | Has Model (used as a proxy) |
| Yuba River | Has Model | Has Model | Has Model |
ws_proxy_summary |>
ggplot() +
facet_wrap(~hab, nrow = 1) +
geom_sf(data = filter(cv_mainstems, habitat != "historical"), color = "lightgray", linetype = "dotted") +
geom_sf(aes(color = label)) +
theme(legend.position = "top") +
guides(color = guide_legend(nrow = 1)) +
scale_color_manual(name = "",
values = c("Has Model" = "#8cc2ca",
"Has Model (used as a proxy)" = "#6388b4",
"Does Not Have Model (uses a proxy)" = "#ef6f6a"))
knitr::knit_exit()